Regresión Lineal Múltiple
TP Integrador
Alumnos
Nicolás Dominutti, Carlos Suarez Gurruchaga, Hernán Telechea
set.seed(101)
COVID = read.csv("COVID.txt", sep="")
dim(COVID)## [1] 139 35summary(COVID)## geoId CntrName muertes casos
## Length:139 Length:139 Min. : 0 Min. : 2.0
## Class :character Class :character 1st Qu.: 12 1st Qu.: 725.5
## Mode :character Mode :character Median : 66 Median : 3333.0
## Mean : 2518 Mean : 37621.0
## 3rd Qu.: 548 3rd Qu.: 20625.0
## Max. :100442 Max. :1699933.0
## muertes.permil casos.permil l10muertes.permil Hombres80
## Min. : 0.000 Min. : 0.949 Min. :0.0000 Min. :0.08737
## 1st Qu.: 2.214 1st Qu.: 70.054 1st Qu.:0.5062 1st Qu.:0.37722
## Median : 6.074 Median : 277.583 Median :0.8497 Median :1.07801
## Mean : 52.735 Mean : 1133.633 Mean :1.0282 Mean :1.59580
## 3rd Qu.: 33.787 3rd Qu.: 1596.139 3rd Qu.:1.5414 3rd Qu.:2.57442
## Max. :819.817 Max. :17596.220 Max. :2.9142 Max. :6.16625
## Mujeres80 Pobla80 Pobla65 PoblaMid
## Min. : 0.1639 Min. :0.1308 Min. : 1.085 Min. :47.42
## 1st Qu.: 0.5983 1st Qu.:0.4885 1st Qu.: 3.284 1st Qu.:60.15
## Median : 1.4421 Median :1.2710 Median : 6.450 Median :64.70
## Mean : 2.6441 Mean :2.1200 Mean : 9.223 Mean :63.58
## 3rd Qu.: 4.5155 3rd Qu.:3.7069 3rd Qu.:14.920 3rd Qu.:66.90
## Max. :10.5564 Max. :8.3613 Max. :27.576 Max. :85.09
## PoblaData PoblaDens Mujeres Urbano
## Min. : 0.001102 Min. : 0.02041 Min. :24.50 Min. : 13.03
## 1st Qu.: 0.047427 1st Qu.: 0.33375 1st Qu.:49.72 1st Qu.: 44.62
## Median : 0.111232 Median : 0.80729 Median :50.40 Median : 62.45
## Mean : 0.485381 Mean : 1.87767 Mean :49.89 Mean : 60.83
## 3rd Qu.: 0.333277 3rd Qu.: 1.46399 3rd Qu.:51.01 3rd Qu.: 78.83
## Max. :13.927300 Max. :79.52998 Max. :54.54 Max. :100.00
## ExpectVida NeontlMort DisMort Lesion
## Min. :49.84 Min. : 0.90 Min. : 7.80 Min. : 2.600
## 1st Qu.:64.57 1st Qu.: 3.40 1st Qu.:14.30 1st Qu.: 5.650
## Median :71.11 Median : 8.70 Median :17.90 Median : 8.700
## Mean :70.06 Mean :12.30 Mean :18.30 Mean : 8.827
## 3rd Qu.:75.14 3rd Qu.:20.45 3rd Qu.:22.35 3rd Qu.:10.700
## Max. :81.70 Max. :42.00 Max. :30.60 Max. :28.400
## EnfNoTrans EnfTrans PBI Tuberculosis
## Min. :26.00 Min. : 1.30 Min. : 0.7568 Min. : 1.0
## 1st Qu.:48.50 1st Qu.: 5.20 1st Qu.: 4.1778 1st Qu.: 13.0
## Median :74.80 Median :11.20 Median : 12.3515 Median : 45.0
## Mean :69.39 Mean :21.78 Mean : 20.0941 Mean :106.5
## 3rd Qu.:88.75 3rd Qu.:38.00 3rd Qu.: 28.3214 3rd Qu.:150.0
## Max. :95.20 Max. :65.30 Max. :123.2139 Max. :611.0
## Diabetes Medicos Camas ImmunSaramp
## Min. : 1.000 Min. :0.01767 Min. : 0.200 Min. :30.00
## 1st Qu.: 5.050 1st Qu.:0.28709 1st Qu.: 1.147 1st Qu.:85.00
## Median : 6.500 Median :1.38303 Median : 2.171 Median :93.00
## Mean : 7.425 Mean :1.65391 Mean : 2.980 Mean :88.21
## 3rd Qu.: 9.200 3rd Qu.:2.77551 3rd Qu.: 3.984 3rd Qu.:97.00
## Max. :22.000 Max. :5.71011 Max. :13.796 Max. :99.00
## TempMarzo HipTen.H HipTen.M HipTen
## Min. :-18.72 Min. : 9.00 Min. :12.40 Min. :10.70
## 1st Qu.: 4.93 1st Qu.:25.85 1st Qu.:25.85 1st Qu.:26.07
## Median : 19.99 Median :32.60 Median :32.00 Median :32.00
## Mean : 15.38 Mean :32.80 Mean :31.98 Mean :32.39
## 3rd Qu.: 25.25 3rd Qu.:40.40 3rd Qu.:38.00 3rd Qu.:39.83
## Max. : 30.63 Max. :54.60 Max. :50.60 Max. :50.60
## BCG BCGf Tiempo
## Min. :0.0000 Length:139 Min. : 0.00
## 1st Qu.:1.0000 Class :character 1st Qu.: 57.00
## Median :1.0000 Mode :character Median : 65.00
## Mean :0.9496 Mean : 63.15
## 3rd Qu.:1.0000 3rd Qu.: 75.00
## Max. :1.0000 Max. :136.00Hay 35 variables analizadas para las 139 filas de datos. Una de las variables es la lista de países. 3 de las variables son cualitativas, mientras que las demás son cuantitativas.
attach(COVID)
par(mfrow=c(1, 2))
boxplot(muertes.permil ~ BCGf,
main = "Muertes por millón de habitantes",
ylab = "Muertes.permil", col=c("thistle","wheat"),
cex = 0.7, xlab = "Estrategia de inmunización",
names = c("selectiva", "todos"))
boxplot(l10muertes.permil ~ BCGf,
main = "Muertes por millón [log]",
ylab = "log10muertes.permil", col=c("thistle","wheat"),
cex = 0.7, xlab = "Estrategia de inmunización",
names = c("selectiva", "todos"))
Ante esta comparación de boxplots, conviene considerar la variables l10muertes.permil porque la relación entre BCGf y muertes.permil no parece ser muy lineal y estaría considerando muchos valores como potencialmente atípicos cuando en realidad no lo son. Al aplicar logaritmo, apenas un dato es considerado extraño, mientras que el resto de los datos parecen estar más claramente distribuidos. Me doy cuenta mucho más fácilmente que tener una estrategia de vacunación contra BCG selectiva podría generar un mayor número de muertes en las poblaciones.
Continuamos el análisis con la variable l10muertes.permil.
#cor(COVID[,-c(1,2,34)])Las relaciones más fuertes entre variables son: Pobla80 con Hombre80 y Mujer80, que tiene mucho sentido dado que la población es una función de los hombres y mujeres mayores a 80 años. El total de mujeres en la población total es mayor a la de los hombres, por eso su correlación es ~0.99 (la de los hombres es ~0.98). En cambio, las muertes no parecen tener relación alguna con otras enfermedades o la densidad poblacional (correlaciones menores al módulo de 0.05).
Planteamos algunas preguntas para identificar relaciones entre variables (muchas de ellas dentro de las más consultadas durante la pandemia):
- ¿Mueren más donde peores condiciones económicas hay?
plot(PBI,l10muertes.permil)
cor(PBI,l10muertes.permil)## [1] 0.5385521Es poco clara la relación. Si bien un alto porcentaje de los países con bajo PBI tienen un bajo total de muertes logarítmicas por millón, es poco claro el incremento (o decrecimiento) de esta variable a mayores niveles de PBI. Es débil esta relación: la correlación es positiva en ~0.53.
- ¿Tiene que ver la edad de los casos en la probabilidad de morir?
plot(PoblaMid,l10muertes.permil)
plot(Pobla65,l10muertes.permil)
plot(Pobla80,l10muertes.permil)
Al comparar 3 tipos de población, identificamos que, cuanto mayor es la edad promedio de la población, más relación lineal hay con respecto al total de muertes. Esto se comprueba con las correlaciones:
cor(PoblaMid,l10muertes.permil)## [1] 0.3623718cor(Pobla65,l10muertes.permil)## [1] 0.6579062cor(Pobla80,l10muertes.permil)## [1] 0.6676371Las 3 correlaciones son positivas pero la mayor se da con Pobla80 (0.6676).
- ¿El COVID se propaga en ambientes más fríos?
plot(TempMarzo,l10muertes.permil)
cor(TempMarzo,l10muertes.permil)## [1] -0.5107478Al parecer, temporaturas más elevadas parecen generar menores chances de mortalidad, quizás porque los contagios son menores y/o el virus resiste menos. De todas formas, no es fuerte la relación (correlación es negativa en ~0.51).
mean(l10muertes.permil[TempMarzo > 20] < 1)## [1] 0.7971014El ~80% de los países cuya temperatura promedio en marzo fue mayor a 20 grados tuvieron menos de 1 muerte logarítmica por mil.
plot(Tiempo,TempMarzo)
mean(Tiempo[TempMarzo > 20] > 60)## [1] 0.8405797Sin embargo, estos datos solo son de marzo de 2020, época en la que la enfermedad aún no estaba tan propagada en países del hemisferio sur (más cálido en esa época). Por ende, en lugar de decir que hubo menos muertes porque el ambiente era más cálido, podríamos inferir que esto fue así porque la enfermedad tardó más en llegar. Con esto, también visualizamos que los países con mayores temperaturas en marzo son también los que más tardaron en registrar contagios, en promedio. De todas maneras, no sabemos si esos contagios tardaron en llegar por la temperatura o por alguna otra causa que desconocemos.
- ¿Hay menos muertes en los países que mayor atención médica tienen?
plot(Camas,l10muertes.permil)
cor(Camas,l10muertes.permil)## [1] 0.3636445plot(Medicos,l10muertes.permil)
cor(Medicos,l10muertes.permil)## [1] 0.6191814df = COVID[ ,c(5:7, 10, 11, 12, 23, 24, 25, 26, 27, 35)]
pairs(df[ ,c("l10muertes.permil", "Pobla80", "Pobla65",
"PBI", "Tuberculosis", "Medicos", "Camas")])
Esta relación parece contradictoria. Si bien no es fuerte, mucho no dice que tener camas ayuda a que la gente no se muera. No es un factor tan determinante. La relación entre médicos y muertes tiene un poco más de sentido. En los mismos lugares en los que mayor número de muertes hay, mayor es la presencia de médicos, quizás por el tamaño poblacional. La correlación es positiva en ~0.62.
- ¿Es más fácil la propagación en ambientes urbanos?
plot(Urbano,l10muertes.permil)
cor(Urbano,l10muertes.permil)## [1] 0.5717382Es bastante claro que ambientes más urbanizados tienen mayor número de muertes. Esto es porque la población total también es mayor.
- ¿Cuál es la incidencia de las muertes por sobre los casos?
COVID$Incidencia = round(COVID$muertes / COVID$casos,2)
hist(COVID$Incidencia)
Vemos que, en la mayoría de los países, la incidencia fue de entre 0 y 0.05. Veamos cuáles son los que quedan hacia la cola derecha:
COVID %>%
filter((COVID$Incidencia>.1)) %>%
select('CntrName')Como punto en común, podemos decir que todos estos países pertenecen al hemisferio norte y, salvo México, Belize y Yemen, +70% de ellos son europeos.
df = COVID %>%
mutate(alto = case_when(Incidencia > .1 ~ 1, Incidencia <= .1 ~ 0))
pairs(df %>%
filter(df$alto==1) %>%
select('Pobla80','muertes.permil',
'casos.permil','l10muertes.permil'),
col = "darkblue")
En este pairplot, podemos ver que, dentro de los países donde la mortalidad del virus sobre los casos totales fue mayor al .1%, existe una posible relación marcada entre la proporción de habitantes con +80 años y las muertes por mil.
aux = df %>% filter(df$alto==1)
plot(y=aux$Pobla80, x=aux$muertes.permil, main='Efecto de poblaciones longevas\n sobre muertes en países con letalidad alta')
abline(h=4.5, col='red')
En este deepdive, podemos ver que se podría segmentar estos países en 2 clusters, teniendo en cuenta el marcado salto que se produce en las muertes por mil cuando la variable Pobla80 supera el 4.4%.
Tomamos otras variables que, intuitivamente, parecieran tener relación directa en muertes o, en algunos casos, en los casos de covid.
¿Existirá alguna relación detectable?
pairs(df %>%
select('ExpectVida','Pobla80','l10muertes.permil'), col = "darkblue")
Vemos que, en países donde la expectativa de vida es mayor, se produjeron más muertes. Esto se podría explicar por la relación directa (logarítmica) que tiene esta variable con la proporción de personas mayores a 80 años presentes en este tipo de países y, cómo vimos en el punto anterior, esta variable se encuentra relacionada con las muertes.
Ahora, intentemos detectar estas relaciones lineales de las que hablamos previamente, calculando el índice de correlación de Pearson:
corrplot(
cor(select_if(df, is.numeric)),
type='lower',
#addCoef.col="black",
number.cex=0.75,
mar=c(0,0,1,0),
title = 'Matriz de correlaciones'
)
lm.covid = lm(l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
Tuberculosis + Camas + TempMarzo + PBI,
data = COVID)
lm.covid_2 = lm(muertes.permil ~ PoblaDens + Pobla80 + Urbano +
Tuberculosis + Camas + TempMarzo + PBI,
data = COVID)
summary(lm.covid)##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08088 -0.29012 -0.03314 0.27929 1.29267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5981046 0.2012630 2.972 0.00352 **
## PoblaDens -0.0111873 0.0063386 -1.765 0.07990 .
## Pobla80 0.2040913 0.0338369 6.032 1.55e-08 ***
## Urbano 0.0082702 0.0025767 3.210 0.00167 **
## Tuberculosis -0.0005815 0.0003431 -1.695 0.09246 .
## Camas -0.1101291 0.0263087 -4.186 5.17e-05 ***
## TempMarzo -0.0143117 0.0057068 -2.508 0.01337 *
## PBI 0.0062466 0.0027181 2.298 0.02314 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4765 on 131 degrees of freedom
## Multiple R-squared: 0.6239, Adjusted R-squared: 0.6038
## F-statistic: 31.04 on 7 and 131 DF, p-value: < 2.2e-16summary(lm.covid_2)##
## Call:
## lm(formula = muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
##
## Residuals:
## Min 1Q Median 3Q Max
## -230.52 -42.58 -8.56 21.49 650.26
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.84325 43.17327 0.251 0.802087
## PoblaDens -1.03060 1.35971 -0.758 0.449840
## Pobla80 40.79237 7.25840 5.620 1.1e-07 ***
## Urbano 0.30373 0.55274 0.550 0.583594
## Tuberculosis 0.05273 0.07359 0.717 0.474932
## Camas -22.21224 5.64354 -3.936 0.000134 ***
## TempMarzo -1.59048 1.22418 -1.299 0.196149
## PBI 1.18989 0.58307 2.041 0.043285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 102.2 on 131 degrees of freedom
## Multiple R-squared: 0.3765, Adjusted R-squared: 0.3432
## F-statistic: 11.3 on 7 and 131 DF, p-value: 3.772e-11Optamos por utilizar como variable target l10muertes.permil, debido a su mayor R2 y cantidad de regresoras significativas para el modelo.
leverage_sorted=sort(hatvalues(lm.covid), decreasing = T)
boxplot(leverage_sorted,
main = "LEVERAGE DE LAS OBSERVACIONES PARA DISTINTOS PAISES",
col = "blue")
leverage_sorted[1:10]## SG QA JP MN CA BY LU BN
## 0.9451793 0.2946915 0.2162310 0.1745735 0.1486524 0.1356584 0.1347664 0.1336935
## GA LS
## 0.1333837 0.1307780Se observa que el modelo le asigna un leverage muy elevado a la observacion que representa a Singapur, seguida por las que representan a Qatar, Japan, Mongolia, Canada, Bielorusia, Luxemburgo, Brunei_Darussalam, Gabon y Lesotho, respectivamente.
plot(lm.covid, 4)
En esta representancion de las distancias de Cook, observamos tambien que las observaciones de Singapur, Japan y Qatar, presentan elevadas distancias de Cook, lo cual reafirma que tienen un leverage elevado, ya que la misma es proporcional a el leverage de la observacion, asi como tambien al residuo de la misma. Por lo que ahora queda pendiente analizar los residuos de estas observaciones para descartar o no que sean posibles outliers, o simplemente observaciones muy influyentes para nuestro modelo.
plot(lm.covid, 3)
Al analizar el aspecto de los residuos del modelo, pareciera haber un grado de HETEROCEDASTICIDAD de los residuos, pero no visualizamos ninguna estructura en los mismos.
cor.test(abs(lm.covid$residuals), lm.covid$fitted.values, method = "spearman")##
## Spearman's rank correlation rho
##
## data: abs(lm.covid$residuals) and lm.covid$fitted.values
## S = 367626, p-value = 0.0355
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1786362Con un pvalue pequeño rechazamos H0, y por tanto las varianzas de los residuos son distintas, hay HETEROCEDASTICIDAD.
plot(lm.covid, 2)
ks.test(lm.covid$residuals, "pnorm",
mean(lm.covid$residuals),
sd(lm.covid$residuals))##
## One-sample Kolmogorov-Smirnov test
##
## data: lm.covid$residuals
## D = 0.048945, p-value = 0.8931
## alternative hypothesis: two-sidedCon un pvalue de 0.89, nuestros residuos tienen una distribucion que se asemeja a una normal.
plot(lm.covid, 5)
Al analizar los residuos en funcion del Leverage de las observaciones, vemos que la observacion correspondiente a Singapur, si bien tiene un leverage elevado como mencionamos en el item anterior, no posee un residuo alto, y por tanto (por el momento) preferimos dejarla en el modelo. Tambien vemos resaltadas las observaciones Qatar y Japan, las cuales tambien poseen un alto valor de Leverage, pero su residuo no parece ser alto y por tanto no serian outliers.
step(lm.covid, direction = "both")## Start: AIC=-198.35
## l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis +
## Camas + TempMarzo + PBI
##
## Df Sum of Sq RSS AIC
## <none> 29.738 -198.34
## - Tuberculosis 1 0.6522 30.390 -197.33
## - PoblaDens 1 0.7071 30.445 -197.08
## - PBI 1 1.1989 30.937 -194.85
## - TempMarzo 1 1.4277 31.166 -193.83
## - Urbano 1 2.3385 32.076 -189.82
## - Camas 1 3.9778 33.716 -182.90
## - Pobla80 1 8.2586 37.996 -166.28##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
##
## Coefficients:
## (Intercept) PoblaDens Pobla80 Urbano Tuberculosis
## 0.5981046 -0.0111873 0.2040913 0.0082702 -0.0005815
## Camas TempMarzo PBI
## -0.1101291 -0.0143117 0.0062466Realizamos un analisis de los residuos del modelo en funcion de las covariables utilizadas a partir del criterio “STEPWISE” y llegamos a la conclusion que el mejor modelo es el que utiliza todas las variables propuestas en el item 3), sin quitar ninguna.
Respecto a las covariables que mas influencia cobran para predecir el target, se encuentran el promedio de poblacion masculina y femenina mayor a 80 años en primer medida, seguido por la cantidad de camas de hospital cada 1000 personas.
ANALISIS DE OUTLIERS POR REGRESION ROBUSTA
Vamos a plantear una regresion robusta para chequear si los betas estimados por OLS dan parecidos, o hay observaciones que nos estan afectando la regresion y no las encontramos en los analisis hechos.
control = lmrobdet.control(bb = 0.5, efficiency = 0.85, family = "bisquare")
modelo_MM_covid = lmrobdetMM(l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
Tuberculosis + Camas + TempMarzo + PBI,
data = COVID, control = control)
summary(modelo_MM_covid)##
## Call:
## lmrobdetMM(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis +
## Camas + TempMarzo + PBI, data = COVID, control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -1.24220 -0.27260 -0.03026 0.31517 1.31445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3849847 0.2222787 1.732 0.08563 .
## PoblaDens -0.0125395 0.0063513 -1.974 0.05045 .
## Pobla80 0.2524079 0.0373247 6.762 4.07e-10 ***
## Urbano 0.0091258 0.0027644 3.301 0.00124 **
## Tuberculosis -0.0004969 0.0003576 -1.390 0.16697
## Camas -0.1163727 0.0286581 -4.061 8.37e-05 ***
## TempMarzo -0.0091945 0.0063197 -1.455 0.14809
## PBI 0.0060534 0.0029396 2.059 0.04145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.4764
## Multiple R-squared: 0.6361, Adjusted R-squared: 0.6166
## Convergence in 25 IRWLS iterations(lm.covid$coefficients)## (Intercept) PoblaDens Pobla80 Urbano Tuberculosis
## 0.5981046383 -0.0111873307 0.2040912643 0.0082702004 -0.0005814871
## Camas TempMarzo PBI
## -0.1101291038 -0.0143117235 0.0062466456(modelo_MM_covid$coefficients)## (Intercept) PoblaDens Pobla80 Urbano Tuberculosis
## 0.3849846701 -0.0125395234 0.2524078877 0.0091257877 -0.0004969057
## Camas TempMarzo PBI
## -0.1163726595 -0.0091945121 0.0060533946Los beta estimados por ambos metodos son muy semejantes, lo que nos da la pauta de que NO tenemos outliers que esten viciando nuestro modelo original de LS, respecto a los coeficientes que no dan estadisiticamente significativos para el metodo robusto, puede deberse a que el mismo por el metodo que emplea para su calculo, sus intervalos son mas largos que los que devuelve OLS y por tanto, puede ocurrir que el 0 caiga dentro del intervalo, dandonos a creer que una variable no tiene informacion util para darnos, cuando no es asi. Por esto vamos a continuar con el metodo de OLS, dado que no tenemos outliers y es el metodo que nos devuelve los intervalos de los betas mas pequeños.
plot(y = modelo_MM_covid$rweights,
x = modelo_MM_covid$residuals, col = "blue",
ylim = c(0, 1), xlab = "RESIDUOS", ylab = "RWEIGHTS",
main = "RWEIGHTS vs RESIDUOS MODELO ROBUSTO")
abline(h = 0, col = "red", lwd = 2)
Reafirmamos que no tenemos observaciones que sean outliers, ya que el metodo MMestimadores, no tiene ninguna observacion que le asigne un peso de 0
df2 = COVID[, c("PoblaDens","Pobla80", "Urbano", "Tuberculosis", "Camas", "TempMarzo", "PBI")]
pairs(df2)
Probamos crear una interaccion entre las variables PBI y Camas
df2[, "PBIcamas"] = df2$Camas * df2$PBI
lm.covid_interaccion = lm(l10muertes.permil ~ PoblaDens + Pobla80 +
Urbano + Tuberculosis + Camas +
TempMarzo + PBI + PBIcamas, data = df2)
summary(lm.covid_interaccion)##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI + PBIcamas, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.07242 -0.29205 -0.03691 0.28372 1.29236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6201428 0.2148048 2.887 0.00456 **
## PoblaDens -0.0111062 0.0063665 -1.744 0.08344 .
## Pobla80 0.2004483 0.0360438 5.561 1.46e-07 ***
## Urbano 0.0084492 0.0026531 3.185 0.00181 **
## Tuberculosis -0.0005942 0.0003468 -1.713 0.08907 .
## Camas -0.1181358 0.0374609 -3.154 0.00200 **
## TempMarzo -0.0147272 0.0058904 -2.500 0.01366 *
## PBI 0.0051970 0.0044247 1.175 0.24233
## PBIcamas 0.0003444 0.0011433 0.301 0.76369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4781 on 130 degrees of freedom
## Multiple R-squared: 0.6241, Adjusted R-squared: 0.601
## F-statistic: 26.98 on 8 and 130 DF, p-value: < 2.2e-16NO observamos una mejoria del R2 respecto al modelo original.
df2[,"PBIpobla80"] = df2$PBI * df2$Pobla80
lm.covid_interaccion_2 = lm(l10muertes.permil ~ PoblaDens + Pobla80 +
Urbano + Tuberculosis + Camas +
TempMarzo + PBI + PBIpobla80, data = df2)
summary(lm.covid_interaccion_2)##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI + PBIpobla80, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08353 -0.28711 -0.05352 0.31835 1.30154
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6286418 0.2019852 3.112 0.002283 **
## PoblaDens -0.0118262 0.0063384 -1.866 0.064322 .
## Pobla80 0.1436949 0.0565796 2.540 0.012272 *
## Urbano 0.0089970 0.0026267 3.425 0.000823 ***
## Tuberculosis -0.0006721 0.0003488 -1.927 0.056176 .
## Camas -0.1048197 0.0265341 -3.950 0.000127 ***
## TempMarzo -0.0139004 0.0056986 -2.439 0.016066 *
## PBI 0.0031349 0.0035807 0.875 0.382926
## PBIpobla80 0.0016874 0.0012690 1.330 0.185937
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4751 on 130 degrees of freedom
## Multiple R-squared: 0.6289, Adjusted R-squared: 0.6061
## F-statistic: 27.54 on 8 and 130 DF, p-value: < 2.2e-16Tampoco observamos un aumento significativo en el R2 al interaccionar PBI con Pobla80.
Continuamos probando interacciones en busca de una mejora del R2.
df2[, "PBItemp"] = df2$TempMarzo * df2$PBI
lm.covid_interaccion_3 = lm(l10muertes.permil ~ PoblaDens + Pobla80 +
Urbano + Tuberculosis + Camas +
TempMarzo + PBI + PBItemp, data = df2)
summary(lm.covid_interaccion_3)##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI + PBItemp, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.03547 -0.28698 -0.06073 0.30877 1.30492
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4306473 0.2233762 1.928 0.056048 .
## PoblaDens -0.0075220 0.0066628 -1.129 0.260995
## Pobla80 0.1820208 0.0360832 5.044 1.5e-06 ***
## Urbano 0.0089768 0.0025934 3.461 0.000728 ***
## Tuberculosis -0.0006461 0.0003429 -1.884 0.061764 .
## Camas -0.0994583 0.0268896 -3.699 0.000319 ***
## TempMarzo -0.0067999 0.0072204 -0.942 0.348057
## PBI 0.0113280 0.0040551 2.793 0.006003 **
## PBItemp -0.0003461 0.0002061 -1.679 0.095517 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4732 on 130 degrees of freedom
## Multiple R-squared: 0.6319, Adjusted R-squared: 0.6092
## F-statistic: 27.89 on 8 and 130 DF, p-value: < 2.2e-16df2[, "PBIurbano"] = df2$Urbano * df2$PBI
lm.covid_interaccion_4 = lm(l10muertes.permil ~ PoblaDens + Pobla80 +
Urbano + Tuberculosis + Camas +
TempMarzo + PBI + PBIurbano, data = df2)
summary(lm.covid_interaccion_4)##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI + PBIurbano, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05951 -0.28139 -0.02248 0.31277 1.27442
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4882083 0.2176972 2.243 0.026616 *
## PoblaDens -0.0086414 0.0066163 -1.306 0.193831
## Pobla80 0.1931678 0.0347706 5.555 1.50e-07 ***
## Urbano 0.0096886 0.0027905 3.472 0.000702 ***
## Tuberculosis -0.0005501 0.0003430 -1.604 0.111162
## Camas -0.1138959 0.0263971 -4.315 3.14e-05 ***
## TempMarzo -0.0140060 0.0056964 -2.459 0.015258 *
## PBI 0.0230117 0.0131379 1.752 0.082208 .
## PBIurbano -0.0001902 0.0001458 -1.304 0.194486
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4752 on 130 degrees of freedom
## Multiple R-squared: 0.6287, Adjusted R-squared: 0.6059
## F-statistic: 27.52 on 8 and 130 DF, p-value: < 2.2e-16df2[, "PBI^2"] = df2$PBI * df2$PBI
lm.covid_interaccion_5 = lm(l10muertes.permil ~ PoblaDens + Pobla80 +
Urbano + Tuberculosis + Camas +
TempMarzo + PBI + PBI^2, data = df2)
summary(lm.covid_interaccion_5)##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI + PBI^2, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08088 -0.29012 -0.03314 0.27929 1.29267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5981046 0.2012630 2.972 0.00352 **
## PoblaDens -0.0111873 0.0063386 -1.765 0.07990 .
## Pobla80 0.2040913 0.0338369 6.032 1.55e-08 ***
## Urbano 0.0082702 0.0025767 3.210 0.00167 **
## Tuberculosis -0.0005815 0.0003431 -1.695 0.09246 .
## Camas -0.1101291 0.0263087 -4.186 5.17e-05 ***
## TempMarzo -0.0143117 0.0057068 -2.508 0.01337 *
## PBI 0.0062466 0.0027181 2.298 0.02314 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4765 on 131 degrees of freedom
## Multiple R-squared: 0.6239, Adjusted R-squared: 0.6038
## F-statistic: 31.04 on 7 and 131 DF, p-value: < 2.2e-16Si bien no conseguimos una mejora muy significativa, la mejor interaccion que conseguimos fue entre PBI y Urbano. Vamos a aplicar nuevamente el criterio de seleccion STEPWISE, para ver si adoptamos esta interaccion o mantenemos nuestro modelo origil.
step(lm.covid_interaccion_4, direction = "both")## Start: AIC=-198.15
## l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis +
## Camas + TempMarzo + PBI + PBIurbano
##
## Df Sum of Sq RSS AIC
## - PBIurbano 1 0.3840 29.738 -198.34
## - PoblaDens 1 0.3852 29.739 -198.34
## <none> 29.354 -198.15
## - Tuberculosis 1 0.5809 29.935 -197.43
## - PBI 1 0.6927 30.047 -196.91
## - TempMarzo 1 1.3650 30.719 -193.83
## - Urbano 1 2.7219 32.076 -187.83
## - Camas 1 4.2036 33.557 -181.55
## - Pobla80 1 6.9689 36.323 -170.54
##
## Step: AIC=-198.35
## l10muertes.permil ~ PoblaDens + Pobla80 + Urbano + Tuberculosis +
## Camas + TempMarzo + PBI
##
## Df Sum of Sq RSS AIC
## <none> 29.738 -198.34
## + PBIurbano 1 0.3840 29.354 -198.15
## - Tuberculosis 1 0.6522 30.390 -197.33
## - PoblaDens 1 0.7071 30.445 -197.08
## - PBI 1 1.1989 30.937 -194.85
## - TempMarzo 1 1.4277 31.166 -193.83
## - Urbano 1 2.3385 32.076 -189.82
## - Camas 1 3.9778 33.716 -182.90
## - Pobla80 1 8.2586 37.996 -166.28##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI, data = df2)
##
## Coefficients:
## (Intercept) PoblaDens Pobla80 Urbano Tuberculosis
## 0.5981046 -0.0111873 0.2040913 0.0082702 -0.0005815
## Camas TempMarzo PBI
## -0.1101291 -0.0143117 0.0062466Nuestro analisis de residuos por STEPWISE considera que es mejor no tomar en cuenta la interaccion, asi que vamos a continuar con nuestro modelo original.
int_prediccion = predict(lm.covid, newdata = COVID, interval = c("predict"))
plot(y = lm.covid$fitted.values,
x = COVID$l10muertes.permil,
col = "darkblue", ylim = c(-2, 5), cex = .8,
xlab = "VALORES OBSERVADOS", ylab = "VALORES PREDICHOS")
# Intervalos de prediccion
points(x = COVID$l10muertes.permil,
y = int_prediccion[, 2], col = "darkred", lwd =3, cex = 0.5)
points(x = COVID$l10muertes.permil,
y = int_prediccion[, 3], col = "darkred", lwd =3, cex = 0.5)
legend("topleft", c("INTERVALO DE PREDICCION", "OBSERVACION"), fill = c("darkred", "darkblue"))
PROBLEMAS DE COLINEALIDAD:
Por construccion sabemos que geoID nos otorga la misma informacion que CntrName.
Por construccion sabemos que casos y casos.permil, estan altamente correlacionadas.
Por construccion sabemos que muertes y muertes.permil, estan altamente correlacionadas.
Por construccion sabemos que Pobla80, esta relacionado altamente con Mujeres80 y Hombres80, ya que es el promedio de los dos anteriores.
Por construccion sabemos tambien que EnfNoTrans y EnfTrans guardan una colinealidad por construccion.
Por construccion sabemos tambien que HipTen esta relacionada altamente con HipTen.H e HipTen.M, ya que es el promedio de los dos anteriores.
Por construccion sabemos tambien que BCG y BCGf tienen correlacion 1.
Sobre estos casos, una posible solucion seria optar por alguna de las variables que estan correlacionadas entre si, eliminando las demas, que no aportan informacion util nueva. Por lo tanto vamos a adoptar quedarnos con CntrName, casos.permil, muertes.permil, Pobla80, EnfTrans, HipTen, BCG, eliminando las demas.
A su vez, tambien vamos a realizar un analisis de multicolinealidad a partir del Factor de inflacion de la varianza (VIF)
df_for_vif = COVID[, -c(1, 2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
lm_for_vif = lm(l10muertes.permil ~., data = df_for_vif)
barplot(vif(lm_for_vif)[vif(lm_for_vif)>5], main = "VIF DE LAS COVARIABLES",
col = "blue", cex.names = 0.65)
abline(h = 5, lwd = 2, col = "red")
Recordemos que el VIF es una medida de inflación de los coeficientes de la variable independiente debido a multicolinealidad entre las variables regresoras. Un VIF que excede 5 requiere investigación de multicolinealidad y un VIF de >10 indica multicolinealidad.
Investigaremos luego los siguientes casos de alto VIF más elevados:
- Pobla80
- Pobla65
- ExpectVida
- EnfTrans
cov_altoVIF = names(vif(lm_for_vif)[vif(lm_for_vif)>5])
data.frame(t(cor(df_for_vif$Pobla80, df_for_vif[, cov_altoVIF],
method = "spearman")))data.frame(t(cor(df_for_vif$Pobla65, df_for_vif[, cov_altoVIF],
method = "spearman")))data.frame(t(cor(df_for_vif$ExpectVida, df_for_vif[, cov_altoVIF],
method = "spearman")))data.frame(t(cor(df_for_vif$EnfTrans, df_for_vif[, cov_altoVIF],
method = "spearman")))- Estos valores pueden deberse a que paises con PBIs mas altos, tienen mejor calidad de vida, con lo que reducen la probabilida de muerte por condiciones nutricionales, ademas poseen mejor infraestructura hospitalaria, con lo que las madres tienen mayores tratamientos disponibles para evitar las enfermedades transmisibles.
- A su vez, paises con mejor infraestructura y mas tratamientos disponibles, aumentan la expectativa de vida de la poblacion, asi como tambien permiten una reduccion de la mortalidad neonatal. Paises con altos PBI, que invierten mas en infraestructura medica, y tienen mayor espectativa de vida, fomentan mas la creacion y atraccion de profesionales en el area de la medicina.
- Deberiamos estudiar como fueron tomados los datos, pero por la tan alta correlacion, podemos llegar a pensar que una persona tiene 80 años, ya tiene mas de 65, con lo cual, la poblacion mayor a 80, esta incluida dentro de la poblacion de mas de 65.
- Si cambiamos el foco a a varible de ExpectVida, vemos que resulta racional que, a mayor expectativa de vida, las variables que indican mortalidad tienden a disminuír, el PBI tiende a aumentar junto con la variable de Médicos cada 1000 habitantes (como vimos antes)
- Por último, cuando miramos las enfermedades transmisibles la historia no parece cambiar y, como distintivo, vemos claramente que esta variable está relacionada con la mortalidad neonatal.
df3 = COVID[, -c(1, 3, 4, 8, 9, 21, 30, 31, 34, 36)]La covariable CntrName, podriamos codificarla como una variable dummie para poder alimentar el modelo, pero optamos por no incluirla en la regresion, ya que realizar ese proceso, aumentaria aun mas la dimensionalidad de nuestro problema que ya de por si es alta, y relacion costo beneficio, no creemos que lo valga.
covariables_l10muertes=df3[, -c(1, 2)]
covariables_muertes_permil=df3[, -c(1, 4)]
modelo_muertes_permil = lm(muertes.permil~. , data = covariables_muertes_permil)
summary(modelo_muertes_permil)##
## Call:
## lm(formula = muertes.permil ~ ., data = covariables_muertes_permil)
##
## Residuals:
## Min 1Q Median 3Q Max
## -225.55 -26.13 -5.17 20.92 354.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -14.32051 607.74146 -0.024 0.98124
## casos.permil 0.02909 0.00578 5.033 1.80e-06 ***
## Pobla80 63.45564 22.45560 2.826 0.00556 **
## Pobla65 -15.27977 7.39576 -2.066 0.04107 *
## PoblaMid -2.86408 3.15996 -0.906 0.36664
## PoblaData 1.12118 4.86447 0.230 0.81813
## PoblaDens -2.08770 1.27482 -1.638 0.10423
## Mujeres 7.12569 4.30409 1.656 0.10054
## Urbano -0.05985 0.53843 -0.111 0.91168
## ExpectVida 2.53633 4.25542 0.596 0.55233
## NeontlMort 1.06751 1.88626 0.566 0.57254
## DisMort 0.05534 3.12626 0.018 0.98591
## Lesion -3.35844 2.71100 -1.239 0.21793
## EnfTrans -1.21433 1.88601 -0.644 0.52095
## PBI 0.07827 0.74049 0.106 0.91600
## Tuberculosis 0.02925 0.08207 0.356 0.72217
## Diabetes -2.95779 2.61345 -1.132 0.26009
## Medicos -10.35879 12.50693 -0.828 0.40925
## Camas -5.09684 5.41314 -0.942 0.34839
## ImmunSaramp 0.18782 0.74213 0.253 0.80066
## TempMarzo 0.85855 1.27358 0.674 0.50158
## HipTen -0.14586 1.23836 -0.118 0.90645
## BCG -209.21052 37.58929 -5.566 1.73e-07 ***
## Tiempo -0.60811 0.50497 -1.204 0.23096
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 81.8 on 115 degrees of freedom
## Multiple R-squared: 0.6494, Adjusted R-squared: 0.5793
## F-statistic: 9.261 on 23 and 115 DF, p-value: < 2.2e-16modelo_l10muertes = lm(l10muertes.permil~. , data = covariables_l10muertes)
summary(modelo_l10muertes)##
## Call:
## lm(formula = l10muertes.permil ~ ., data = covariables_l10muertes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.52890 -0.22590 0.00776 0.25171 0.98278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.015e+00 3.230e+00 -1.243 0.21634
## casos.permil 1.670e-04 3.071e-05 5.436 3.1e-07 ***
## Pobla80 6.175e-02 1.193e-01 0.517 0.60584
## Pobla65 2.951e-02 3.930e-02 0.751 0.45423
## PoblaMid 1.632e-02 1.679e-02 0.972 0.33323
## PoblaData 4.407e-03 2.585e-02 0.170 0.86492
## PoblaDens -1.879e-02 6.774e-03 -2.773 0.00647 **
## Mujeres 5.051e-02 2.287e-02 2.208 0.02922 *
## Urbano 6.648e-03 2.861e-03 2.323 0.02192 *
## ExpectVida 2.247e-02 2.261e-02 0.994 0.32244
## NeontlMort 5.602e-03 1.002e-02 0.559 0.57736
## DisMort 5.238e-04 1.661e-02 0.032 0.97490
## Lesion -2.224e-03 1.441e-02 -0.154 0.87761
## EnfTrans 7.238e-03 1.002e-02 0.722 0.47164
## PBI -2.326e-03 3.935e-03 -0.591 0.55567
## Tuberculosis -7.972e-04 4.361e-04 -1.828 0.07016 .
## Diabetes 1.209e-02 1.389e-02 0.871 0.38575
## Medicos 2.294e-02 6.646e-02 0.345 0.73062
## Camas -6.068e-02 2.877e-02 -2.109 0.03708 *
## ImmunSaramp -7.843e-03 3.944e-03 -1.989 0.04911 *
## TempMarzo -7.626e-03 6.768e-03 -1.127 0.26217
## HipTen -1.524e-03 6.581e-03 -0.232 0.81731
## BCG -2.400e-01 1.998e-01 -1.201 0.23204
## Tiempo 4.089e-05 2.683e-03 0.015 0.98787
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4347 on 115 degrees of freedom
## Multiple R-squared: 0.7252, Adjusted R-squared: 0.6702
## F-statistic: 13.19 on 23 and 115 DF, p-value: < 2.2e-16Observamos que utilizando la transformacion logaritmica del target obtenemos el R2 mas elevado, asi como tambien mas estimadores estadisticamente distintos de cero. Aun asi, observamos que debido a la gran varianza que tienen nuestros estimadores, hay muchos coeficientes que no logran ser significativos, esto puede deberse a los problemas de colinealidad que mencionamos previamente y aun no eliminamos.
Mejor modelo hasta el momento: modelo_l10muertes.
Vamos a aplicar el metodo de seleccion de variables stepwise, utilizando como metrica a minimizar el criterio de AKAIKE, para intentar disminuir los problemas de colinealidad y con ello mejorar nuestro modelo.
wise_modelo_l10muertes = step(modelo_l10muertes, direction = c("both"),
trace = F)
summary(wise_modelo_l10muertes)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla80 + PoblaDens +
## Mujeres + Urbano + ExpectVida + Tuberculosis + Camas + ImmunSaramp +
## BCG, data = covariables_l10muertes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.57911 -0.26035 -0.00065 0.26120 0.91573
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.292e+00 1.095e+00 -2.093 0.03831 *
## casos.permil 1.604e-04 2.546e-05 6.299 4.42e-09 ***
## Pobla80 1.373e-01 3.729e-02 3.682 0.00034 ***
## PoblaDens -1.811e-02 5.472e-03 -3.310 0.00121 **
## Mujeres 4.485e-02 1.522e-02 2.946 0.00383 **
## Urbano 5.304e-03 2.370e-03 2.238 0.02695 *
## ExpectVida 2.121e-02 1.056e-02 2.009 0.04666 *
## Tuberculosis -5.820e-04 3.573e-04 -1.629 0.10585
## Camas -4.556e-02 2.185e-02 -2.086 0.03900 *
## ImmunSaramp -7.726e-03 3.489e-03 -2.214 0.02858 *
## BCG -3.009e-01 1.791e-01 -1.680 0.09544 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4219 on 128 degrees of freedom
## Multiple R-squared: 0.7118, Adjusted R-squared: 0.6893
## F-statistic: 31.61 on 10 and 128 DF, p-value: < 2.2e-16barplot(vif(wise_modelo_l10muertes), main = "VIF DE LAS COVARIABLES",
col = "blue", cex.names = 0.65,
ylim=c(0,6), las=2)
abline(h = 5, lwd = 2, col = "red")
Disminuimos la cantidad de covariables de nuestro modelo de 23 a 10, el R2 aumento un poco, el RSE disminuyo tambien y, al parecer, eliminamos el problema de multicolinealidad. A su vez vemos tambien que el numero de betas estimados que son significativamente distintos de cero aumento considerablemente, con lo que tenemos un mejor modelo, con menor varianza en sus coeficientes.
wise_modelo_muertes_permil = step(modelo_muertes_permil, direction = c("both"),
trace = F)
summary(wise_modelo_muertes_permil)##
## Call:
## lm(formula = muertes.permil ~ casos.permil + Pobla80 + Pobla65 +
## PoblaDens + Mujeres + Camas + BCG + Tiempo, data = covariables_muertes_permil)
##
## Residuals:
## Min 1Q Median 3Q Max
## -231.74 -20.62 -1.09 19.98 372.17
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.192e+02 1.388e+02 -1.579 0.116726
## casos.permil 2.801e-02 4.673e-03 5.994 1.89e-08 ***
## Pobla80 6.426e+01 1.776e+01 3.618 0.000423 ***
## Pobla65 -1.414e+01 5.495e+00 -2.574 0.011176 *
## PoblaDens -1.809e+00 1.032e+00 -1.753 0.081909 .
## Mujeres 1.026e+01 2.761e+00 3.717 0.000298 ***
## Camas -9.202e+00 4.020e+00 -2.289 0.023698 *
## BCG -2.091e+02 3.357e+01 -6.227 6.07e-09 ***
## Tiempo -7.640e-01 3.846e-01 -1.986 0.049081 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 78.64 on 130 degrees of freedom
## Multiple R-squared: 0.6337, Adjusted R-squared: 0.6111
## F-statistic: 28.11 on 8 and 130 DF, p-value: < 2.2e-16Observamos nuevamente que conviene trabajar con el target convertido, pero destacamos la mejora en el R2 y RSE respecto al modelo antes de aplicar el criterio de seleccion.
Vamos a probar utilizar un enfoque backward a ver si llegamos a un
resultado distinto
back_modelo_l10muertes = step(modelo_l10muertes, direction = c("back"),
trace = F)
summary(back_modelo_l10muertes)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla80 + PoblaDens +
## Mujeres + Urbano + ExpectVida + Tuberculosis + Camas + ImmunSaramp +
## BCG, data = covariables_l10muertes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.57911 -0.26035 -0.00065 0.26120 0.91573
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.292e+00 1.095e+00 -2.093 0.03831 *
## casos.permil 1.604e-04 2.546e-05 6.299 4.42e-09 ***
## Pobla80 1.373e-01 3.729e-02 3.682 0.00034 ***
## PoblaDens -1.811e-02 5.472e-03 -3.310 0.00121 **
## Mujeres 4.485e-02 1.522e-02 2.946 0.00383 **
## Urbano 5.304e-03 2.370e-03 2.238 0.02695 *
## ExpectVida 2.121e-02 1.056e-02 2.009 0.04666 *
## Tuberculosis -5.820e-04 3.573e-04 -1.629 0.10585
## Camas -4.556e-02 2.185e-02 -2.086 0.03900 *
## ImmunSaramp -7.726e-03 3.489e-03 -2.214 0.02858 *
## BCG -3.009e-01 1.791e-01 -1.680 0.09544 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4219 on 128 degrees of freedom
## Multiple R-squared: 0.7118, Adjusted R-squared: 0.6893
## F-statistic: 31.61 on 10 and 128 DF, p-value: < 2.2e-16Vemos que el R2 y el RSE coinciden con el metodo stepwise, asi como los betas estimados, con lo que concluimos que el mejor modelo que podemos conseguir a este punto es wise_modelo_l10muertes.
METODOS DE REGULARIZACION
Vamos a probar ahora reducir la dimensionalidad de nuestro problema a partir de metodos de regularizacion.
xx = model.matrix(l10muertes.permil ~., data = covariables_l10muertes)[, -1]
lassocv_modelo_l10muertes <- cv.glmnet(xx, y = l10muertes.permil,
alpha = 1,
nfolds = 5,
family = "gaussian",
type.measure = "mse")Obtenemos los mejores valores de lambda para aplicar regularizacion lasso, usando como metrica el error cuadratico medio.
lassocv_modelo_l10muertes$lambda.min## [1] 0.003990394lassocv_modelo_l10muertes$lambda.1se## [1] 0.3470641Ahora aplicamos la regularizacion lasso con el lambda obtenido para ver que covariables sobreviven.
lasso_modelo_l10muertes_min <- glmnet(xx, y = l10muertes.permil,
alpha = 1,
lambda = lassocv_modelo_l10muertes$lambda.min)
lasso_modelo_l10muertes_1sd <- glmnet(xx, y = l10muertes.permil,
alpha = 1,
lambda = lassocv_modelo_l10muertes$lambda.1se)
cbind(coef(lasso_modelo_l10muertes_min),
coef(lasso_modelo_l10muertes_1sd))## 24 x 2 sparse Matrix of class "dgCMatrix"
## s0 s0
## (Intercept) -1.7245118950 2.506694e-01
## casos.permil 0.0001601537 4.219719e-06
## Pobla80 0.0701753098 5.355739e-02
## Pobla65 0.0232764713 .
## PoblaMid 0.0040921465 .
## PoblaData 0.0043040149 .
## PoblaDens -0.0165366840 .
## Mujeres 0.0386786518 .
## Urbano 0.0059785151 .
## ExpectVida 0.0125988089 9.408924e-03
## NeontlMort 0.0040718750 .
## DisMort -0.0039817202 .
## Lesion -0.0053461327 .
## EnfTrans . .
## PBI -0.0009665838 .
## Tuberculosis -0.0006294161 .
## Diabetes 0.0083184355 .
## Medicos . .
## Camas -0.0533819815 .
## ImmunSaramp -0.0067932573 .
## TempMarzo -0.0057097097 .
## HipTen -0.0005048359 .
## BCG -0.2552166366 .
## Tiempo . .Notamos que usando la regularizacion con lambda.min, no logramos reducir demasiado la dimensionalidad del problema, por lo que vamos a adoptar directamente el criterio de seleccion por lambda.1se
CRITERIO DE LAMBDA MIN + 1se
Sobreviven las covariables casos.permil, Pobla80, ExpectVida.
Si bien el metodo nos quita la variable BCG, por el tipo de problema que se esta analizando, consideramos que es importante, con lo que vamos a conservala de todos modos.
ls_covid_lasso_final_1se = lm(l10muertes.permil ~ casos.permil + Pobla80 + ExpectVida +
BCG, data = covariables_l10muertes)
summary(ls_covid_lasso_final_1se)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla80 + ExpectVida +
## BCG, data = covariables_l10muertes)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.78427 -0.31310 0.04202 0.29393 1.17570
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.102e-01 5.543e-01 -0.379 0.7051
## casos.permil 1.222e-04 2.194e-05 5.571 1.34e-07 ***
## Pobla80 1.618e-01 2.751e-02 5.881 3.08e-08 ***
## ExpectVida 1.750e-02 7.897e-03 2.216 0.0284 *
## BCG -4.939e-01 1.912e-01 -2.584 0.0108 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4634 on 134 degrees of freedom
## Multiple R-squared: 0.6361, Adjusted R-squared: 0.6253
## F-statistic: 58.56 on 4 and 134 DF, p-value: < 2.2e-16Notamos que al utilizar el metodo de regularizacion lasso, nuestro R2 NO supera al obtenido a partir del criterio stepwise. Por lo tanto vamos a conservar el modelo obtenido con stepwise, como el definitivo.
Entonces por el momento nuestro mejor modelo es el wise_modelo_l10muertes, vamos a analizar sus residuos para entender la performance
plot(wise_modelo_l10muertes)
Vemos la clara presencia de un outlier (alto leverage y alto residuo) que está dañando la regresión. El país que causa el problema es Qatar y la razón por la que es determinado como outlier es que tiene valores extremos de casos.permil, sin embargo, es uno de los países que menos muertes por mil tuvo. Tal diferencia nos llevaría a pensar que pudo haberse dado un error tanto en la recolección de los datos de casos o de muertes. También cabe destacar que es el país con mayor PBI (aunque este número pareciera ser más racional).
par(mfrow=c(1, 3))
aux = COVID[COVID$casos.permil>3000,]
barplot(aux$casos.permil, names.arg=aux$geoId, las=2, main='Casos por mil por país')
barplot(aux$muertes.permil, names.arg=aux$geoId, las=2, main='Muertes por mil por país')
barplot(aux$PBI, names.arg=aux$geoId, las=2, main='PBI por país')
Tomando este punto de evaluación como disparador, vamos a utilizar distintas técnicas para intentar mejorar la performance de nuestro mejor modelo hasta el momento.
¿Casos atípicos? Proceso de optimización del modelo
Prueba 1: drop de Qatar
Probemos ajustando una regresión sin este valor atípico, la misma debiera performar mucho mejor.
df4 = COVID[, -c(2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
df4 = df4[df4$geoId!='QA',]
modelo_l10muertes_V2 = lm(l10muertes.permil~. , data = df4[,-1])
summary(modelo_l10muertes_V2)##
## Call:
## lm(formula = l10muertes.permil ~ ., data = df4[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73444 -0.21277 -0.01593 0.20441 0.82161
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.770e-01 2.670e+00 -0.291 0.771617
## casos.permil 3.189e-04 3.199e-05 9.969 < 2e-16 ***
## Pobla80 -3.884e-02 9.831e-02 -0.395 0.693488
## Pobla65 6.583e-02 3.244e-02 2.030 0.044728 *
## PoblaMid -1.730e-03 1.391e-02 -0.124 0.901258
## PoblaData 7.795e-03 2.111e-02 0.369 0.712627
## PoblaDens -2.078e-02 5.537e-03 -3.752 0.000277 ***
## Mujeres 2.887e-05 1.980e-02 0.001 0.998839
## Urbano 3.032e-03 2.383e-03 1.272 0.205969
## ExpectVida 1.906e-02 1.847e-02 1.032 0.304132
## NeontlMort 6.926e-03 8.185e-03 0.846 0.399253
## DisMort -1.412e-04 1.356e-02 -0.010 0.991714
## Lesion 1.494e-02 1.197e-02 1.248 0.214567
## EnfTrans 1.734e-04 8.235e-03 0.021 0.983238
## PBI -5.292e-03 3.236e-03 -1.635 0.104735
## Tuberculosis -2.979e-04 3.620e-04 -0.823 0.412187
## Diabetes 1.592e-02 1.135e-02 1.403 0.163449
## Medicos 4.153e-02 5.432e-02 0.765 0.446066
## Camas -4.001e-02 2.364e-02 -1.692 0.093315 .
## ImmunSaramp -6.751e-03 3.223e-03 -2.095 0.038417 *
## TempMarzo -1.475e-03 5.584e-03 -0.264 0.792136
## HipTen 3.105e-03 5.407e-03 0.574 0.566871
## BCG -1.453e-01 1.636e-01 -0.889 0.376128
## Tiempo 1.346e-03 2.197e-03 0.613 0.541383
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3549 on 114 degrees of freedom
## Multiple R-squared: 0.8184, Adjusted R-squared: 0.7818
## F-statistic: 22.34 on 23 and 114 DF, p-value: < 2.2e-16Efectivamente logramos que la regresión aumentara su R2ajustado en +10%, pero vemos que solo 4 de los beta estimados rechazan el test-t, lo que sería un indicio probable de la presencia de multicolinealidad, vamos a correr un proceso de selección de variables sobre este modelo para intentar solucionar dicho problema.
wise_modelo_l10muertes_V2 = step(modelo_l10muertes_V2, direction = c("both"),
trace = F)
summary(wise_modelo_l10muertes_V2)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaDens +
## Urbano + PBI + Tuberculosis + Diabetes + Camas + ImmunSaramp,
## data = df4[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81124 -0.19411 0.01204 0.21291 0.80697
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6287094 0.2507573 2.507 0.0134 *
## casos.permil 0.0003240 0.0000287 11.288 < 2e-16 ***
## Pobla65 0.0589951 0.0074604 7.908 1.06e-12 ***
## PoblaDens -0.0232308 0.0045525 -5.103 1.18e-06 ***
## Urbano 0.0042204 0.0018903 2.233 0.0273 *
## PBI -0.0043099 0.0025914 -1.663 0.0987 .
## Tuberculosis -0.0005758 0.0002664 -2.162 0.0325 *
## Diabetes 0.0156212 0.0087189 1.792 0.0756 .
## Camas -0.0450035 0.0179441 -2.508 0.0134 *
## ImmunSaramp -0.0059518 0.0026293 -2.264 0.0253 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3431 on 128 degrees of freedom
## Multiple R-squared: 0.8094, Adjusted R-squared: 0.796
## F-statistic: 60.4 on 9 and 128 DF, p-value: < 2.2e-16Vemos como al seleccionar las variables correctas pareciera que el modelo ya no tiene multicolinealidades y mejora su performance, pero nos surge una nueva duda ¿existen nuevos outliers que no hayamos detectado hasta el momento?
Prueba 2: nuevos outliers? Drop de Singapur
Como disparador en la respuesta de esta pregunta vamos inspeccionar los residuos del último modelo propuesto
plot(wise_modelo_l10muertes_V2)
Encontramos un nuevo punto con alto leverage y, pese a tener un relativamente bajo residuo estandarizado absoluta, tiene una alta distancia de Cook, por lo que investigaremos más en detalle este punto, que es Singapur (ya visto como punto extremo en el punto 3c del presente trabajo)
aux = COVID[COVID$PoblaDens>3,]
barplot(aux$PoblaDens, names.arg=aux$geoId, las=2, main='Densidad poblacional por país')
Vemos que Singapur es país con mayor densidad poblacional del set de datos, llegando estar a +79 (siendo la segunda mayor de apenas más de 12). Puede ser que, debido a este punto, nuestro modelo no pueda lograr ajustar bien a dicho país (teniendo en cuenta que el beta estimado para la variable PoblaDens es de -0.023). Vamos a probar sacando este valor, volviendo a ajustar y corriendo un proceso de selección de variables para lidiar con la multicolinealidad.
df5 = df4[df4$geoId!='SG',]
modelo_l10muertes_V3 = lm(l10muertes.permil~. , data = df5[,-1])
summary(modelo_l10muertes_V3)##
## Call:
## lm(formula = l10muertes.permil ~ ., data = df5[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73187 -0.21917 -0.01577 0.20766 0.82485
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.408e-01 2.682e+00 -0.276 0.7829
## casos.permil 3.191e-04 3.211e-05 9.939 <2e-16 ***
## Pobla80 -4.196e-02 9.895e-02 -0.424 0.6723
## Pobla65 6.693e-02 3.266e-02 2.049 0.0428 *
## PoblaMid -1.535e-03 1.397e-02 -0.110 0.9127
## PoblaData 6.377e-03 2.146e-02 0.297 0.7668
## PoblaDens -1.206e-02 2.159e-02 -0.559 0.5775
## Mujeres 3.297e-04 1.989e-02 0.017 0.9868
## Urbano 3.270e-03 2.459e-03 1.330 0.1863
## ExpectVida 1.792e-02 1.874e-02 0.956 0.3409
## NeontlMort 6.493e-03 8.280e-03 0.784 0.4346
## DisMort -3.072e-05 1.362e-02 -0.002 0.9982
## Lesion 1.520e-02 1.203e-02 1.263 0.2093
## EnfTrans 3.382e-04 8.274e-03 0.041 0.9675
## PBI -5.204e-03 3.255e-03 -1.599 0.1126
## Tuberculosis -3.092e-04 3.643e-04 -0.849 0.3978
## Diabetes 1.527e-02 1.150e-02 1.328 0.1868
## Medicos 4.259e-02 5.457e-02 0.780 0.4368
## Camas -4.273e-02 2.460e-02 -1.737 0.0852 .
## ImmunSaramp -6.835e-03 3.241e-03 -2.109 0.0372 *
## TempMarzo -2.011e-03 5.749e-03 -0.350 0.7271
## HipTen 2.999e-03 5.432e-03 0.552 0.5820
## BCG -1.232e-01 1.724e-01 -0.715 0.4762
## Tiempo 1.257e-03 2.216e-03 0.567 0.5715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3562 on 113 degrees of freedom
## Multiple R-squared: 0.8184, Adjusted R-squared: 0.7815
## F-statistic: 22.15 on 23 and 113 DF, p-value: < 2.2e-16wise_modelo_muertes_permil_V3 = step(modelo_l10muertes_V3, direction = c("both"),
trace = F)
summary(wise_modelo_muertes_permil_V3)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + Urbano +
## PBI + Tuberculosis + Diabetes + Camas + ImmunSaramp, data = df5[,
## -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81355 -0.20138 0.01161 0.21331 0.81262
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.239e-01 2.504e-01 2.492 0.0140 *
## casos.permil 3.229e-04 2.864e-05 11.276 < 2e-16 ***
## Pobla65 5.849e-02 7.456e-03 7.844 1.49e-12 ***
## Urbano 4.662e-03 1.884e-03 2.475 0.0146 *
## PBI -4.229e-03 2.590e-03 -1.633 0.1050
## Tuberculosis -5.764e-04 2.660e-04 -2.167 0.0321 *
## Diabetes 1.328e-02 8.735e-03 1.521 0.1308
## Camas -4.617e-02 1.793e-02 -2.575 0.0112 *
## ImmunSaramp -6.255e-03 2.624e-03 -2.384 0.0186 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3426 on 128 degrees of freedom
## Multiple R-squared: 0.8097, Adjusted R-squared: 0.7978
## F-statistic: 68.09 on 8 and 128 DF, p-value: < 2.2e-16Vemos que nuestro r2-ajustado mejora un poco y la gran mayoría de las variables resultan significativas
plot(wise_modelo_muertes_permil_V3)
Al ver los residuos podemos decir que ya no se ven puntos atípicos a simple vista, al parecer hemos logrado un modelo consistente y que performa bien (recordemos el salto de un modelo de 0.61 de r2-ajustado a uno de 0.7978). Pero no vamos a detenernos con lo obtenido, tenemos aun 2 approachs por probar, siendo el primero de ellos una regresión robusta
Prueba 3: regresión robusta
Para esta prueba vamos a retrotraernos al modelo antes de dropear cualquier instancia. El plan es simple, ajustar una regresión robusta y luego analizar los residuos de esta, para detectar los outliers.
control = lmrobdet.control(bb = 0.5, efficiency = 0.85, family = "bisquare")
df6 = COVID[, -c(2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
modelo_MM_covid_V2 = lmrobdetMM(l10muertes.permil ~ .,
data = df6[,-1], control = control)
summary(modelo_MM_covid_V2)##
## Call:
## lmrobdetMM(formula = l10muertes.permil ~ ., data = df6[, -1], control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -4.9410 -0.1828 -0.0245 0.1554 0.8794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.748e+00 2.800e+00 0.624 0.533674
## casos.permil 3.253e-04 3.597e-05 9.041 4.44e-15 ***
## Pobla80 -2.054e-02 1.034e-01 -0.199 0.842840
## Pobla65 6.328e-02 3.334e-02 1.898 0.060190 .
## PoblaMid -9.907e-03 1.453e-02 -0.682 0.496568
## PoblaData 3.064e-01 9.013e-02 3.400 0.000927 ***
## PoblaDens -1.981e-02 5.600e-03 -3.537 0.000584 ***
## Mujeres -1.298e-02 2.074e-02 -0.626 0.532507
## Urbano 4.829e-03 2.546e-03 1.897 0.060352 .
## ExpectVida 1.283e-04 2.002e-02 0.006 0.994897
## NeontlMort 2.653e-03 8.600e-03 0.308 0.758299
## DisMort -1.046e-02 1.474e-02 -0.710 0.479353
## Lesion 1.047e-02 1.243e-02 0.842 0.401434
## EnfTrans -1.180e-03 8.424e-03 -0.140 0.888823
## PBI -4.303e-03 3.634e-03 -1.184 0.238827
## Tuberculosis -7.273e-04 3.858e-04 -1.885 0.061924 .
## Diabetes 1.807e-02 1.182e-02 1.529 0.129119
## Medicos -3.670e-02 6.256e-02 -0.587 0.558564
## Camas -1.134e-02 2.785e-02 -0.407 0.684717
## ImmunSaramp -5.581e-03 3.375e-03 -1.653 0.100971
## TempMarzo -5.852e-03 5.937e-03 -0.986 0.326348
## HipTen 3.743e-03 5.562e-03 0.673 0.502362
## BCG -1.680e-01 1.718e-01 -0.978 0.330183
## Tiempo 2.712e-03 2.270e-03 1.195 0.234723
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.3624
## Multiple R-squared: 0.7963, Adjusted R-squared: 0.7555
## Convergence in 36 IRWLS iterationsplot(modelo_MM_covid_V2$residuals,modelo_MM_covid_V2$rweights,
main = "REGRESION ROBUSTA 1", xlab = "RESIDUOS",
ylab = "RWEIGHTS",col = "blue")
Vemos que con este método no vemos 2 valores atípicos como en el anterior, sino que vemos 5 de ellos. Desde nuestro punto de vista resultaría muy “costoso” dropear todos estos valores, por lo que no quisieramos recurrir a esta opción. De todas formas vamos a correr un modelo quitando estos datos para analizar si esto mejora la performance del mismo y, de hacerlo, dimensionar en cuanto lo hace.
df6 = df6[(df6$geoId!='US')&(df6$geoId!='CN')&(df6$geoId!='IN')&(df6$geoId!='QA')&(df6$geoId!='JP'),]# DROPEAMOS LOS 5 OUTLIERS QUE DETECTO MM
modelo_l10muertes_V3 = lm(l10muertes.permil~. , data = df6[,-1])
summary(modelo_l10muertes_V3)##
## Call:
## lm(formula = l10muertes.permil ~ ., data = df6[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69790 -0.18409 -0.00631 0.16924 0.83687
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.403e+00 2.606e+00 0.538 0.591442
## casos.permil 3.127e-04 3.124e-05 10.012 < 2e-16 ***
## Pobla80 -3.494e-02 9.670e-02 -0.361 0.718563
## Pobla65 6.813e-02 3.126e-02 2.179 0.031443 *
## PoblaMid -8.959e-03 1.354e-02 -0.662 0.509491
## PoblaData 2.931e-01 8.317e-02 3.524 0.000620 ***
## PoblaDens -1.919e-02 5.326e-03 -3.603 0.000475 ***
## Mujeres -1.391e-02 1.944e-02 -0.715 0.475819
## Urbano 3.608e-03 2.314e-03 1.559 0.121912
## ExpectVida 3.455e-03 1.845e-02 0.187 0.851793
## NeontlMort 4.801e-03 7.891e-03 0.608 0.544176
## DisMort -7.163e-03 1.315e-02 -0.545 0.587176
## Lesion 1.598e-02 1.145e-02 1.396 0.165653
## EnfTrans -2.993e-03 7.875e-03 -0.380 0.704639
## PBI -4.247e-03 3.112e-03 -1.365 0.175149
## Tuberculosis -4.938e-04 3.520e-04 -1.403 0.163498
## Diabetes 1.302e-02 1.099e-02 1.184 0.238944
## Medicos -1.736e-03 5.629e-02 -0.031 0.975455
## Camas -2.573e-02 2.531e-02 -1.017 0.311615
## ImmunSaramp -4.772e-03 3.106e-03 -1.536 0.127344
## TempMarzo -3.039e-03 5.436e-03 -0.559 0.577186
## HipTen 5.621e-03 5.219e-03 1.077 0.283843
## BCG -2.087e-01 1.631e-01 -1.280 0.203301
## Tiempo 1.904e-03 2.129e-03 0.894 0.373299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3376 on 110 degrees of freedom
## Multiple R-squared: 0.8364, Adjusted R-squared: 0.8021
## F-statistic: 24.44 on 23 and 110 DF, p-value: < 2.2e-16wise_modelo_muertes_permil_V4 = step(modelo_l10muertes_V3, direction = c("both"),
trace = F)
formula_elegida = wise_modelo_muertes_permil_V4$call$formula #Nos guardamos las variables seleccionadas para utilizarlas luego
summary(wise_modelo_muertes_permil_V4)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaData +
## PoblaDens + Urbano + Tuberculosis + Camas + ImmunSaramp +
## BCG, data = df6[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.67764 -0.17239 0.00686 0.16236 0.87887
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.795e-01 2.814e-01 3.125 0.002215 **
## casos.permil 2.902e-04 2.516e-05 11.532 < 2e-16 ***
## Pobla65 5.165e-02 7.238e-03 7.136 7.10e-11 ***
## PoblaData 2.718e-01 7.021e-02 3.871 0.000174 ***
## PoblaDens -2.484e-02 4.317e-03 -5.754 6.43e-08 ***
## Urbano 3.467e-03 1.686e-03 2.056 0.041857 *
## Tuberculosis -7.638e-04 2.510e-04 -3.043 0.002862 **
## Camas -3.255e-02 1.821e-02 -1.787 0.076329 .
## ImmunSaramp -5.360e-03 2.492e-03 -2.150 0.033461 *
## BCG -2.128e-01 1.486e-01 -1.432 0.154605
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3298 on 124 degrees of freedom
## Multiple R-squared: 0.824, Adjusted R-squared: 0.8112
## F-statistic: 64.48 on 9 and 124 DF, p-value: < 2.2e-16Vemos que el modelo ajustado por OLS mejora en un 1.6% (0.7978 vs 0.8112). Vamos a volver a ajustar una regresión robusta para inspeccionar nuevamente el gráfico de residuos que disparó este approach, para ello, vamos a usar las mismas variables que surgieron del último proceso de selección de variables.
modelo_MM_covid_V2 = lmrobdetMM(formula_elegida,
data = select_if(df6,is.numeric), control = control)
summary(modelo_MM_covid_V2)##
## Call:
## lmrobdetMM(formula = formula_elegida, data = select_if(df6, is.numeric),
## control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -0.73780 -0.16801 0.01109 0.18404 0.86496
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.229e-01 3.031e-01 3.045 0.002843 **
## casos.permil 3.058e-04 2.967e-05 10.307 < 2e-16 ***
## Pobla65 5.021e-02 8.032e-03 6.251 6.0e-09 ***
## PoblaData 2.806e-01 7.805e-02 3.595 0.000466 ***
## PoblaDens -2.533e-02 4.468e-03 -5.670 9.5e-08 ***
## Urbano 3.517e-03 1.873e-03 1.877 0.062859 .
## Tuberculosis -8.000e-04 2.715e-04 -2.947 0.003834 **
## Camas -2.320e-02 2.092e-02 -1.109 0.269644
## ImmunSaramp -6.583e-03 2.701e-03 -2.437 0.016209 *
## BCG -1.824e-01 1.565e-01 -1.166 0.246015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.3118
## Multiple R-squared: 0.8243, Adjusted R-squared: 0.8116
## Convergence in 32 IRWLS iterationsplot(modelo_MM_covid_V2$residuals,modelo_MM_covid_V2$rweights,
main = "REGRESION ROBUSTA 2", xlab = "RESIDUOS",
ylab = "RWEIGHTS",col = "blue")
Vemos que tanto los beta no varían mucho contra el modelo de OLS y en el gráfico de rweights vs residuals tenemos una distribución que no indica presencia de outliers. Este termina siendo nuestro mejor modelo por el momento pero, como desventaja, debemos comprender que dropeamos 5 países (Canadá, Qatar, Estados Unidos, India y Japón) podría tener su costo al momento de intentar extrapolar algún tipo de insight que surja del presente trabajo. De todas formas, aun nos queda un approach por probar, el mismo consiste en buscar alguna transformación de la variable target o las covariables regresoras que nos permita subsanar el problema de outliers e intentar mejorar la performance del modelo.
Prueba 4: transformaciones
Como vimos en puntos anteriores una variable que presenta valores extremos es la de casos.permil, por ello resulta lógico probar sobre ella una transformación del tipo logarítmica pero, para corroborar nuestra hipótesis, vamos a correr un test de Box-Tidwell sobre ella para ver qué tipo de transformación nos sugiere
df7 = COVID[, -c(2, 3, 4, 5, 8, 9, 21, 30, 31, 34, 36)]
car::boxTidwell(l10muertes.permil~casos.permil, data=df7[,-1])## MLE of lambda Score Statistic (z) Pr(>|z|)
## 0.125 -14.734 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## iterations = 10Efectivamente, el test Box-Tidwell nos sugiere hacer una transformación logarítmica (dado que un lambda=0 cae dentro de sus intervalos de confianza)
df7$casos.permil = log(df7$casos.permil)
modelo_l10muertes_outliers = lm(l10muertes.permil~. , data = df7[,-1])
summary(modelo_l10muertes_outliers)##
## Call:
## lm(formula = l10muertes.permil ~ ., data = df7[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.78821 -0.19981 0.01797 0.15267 0.71851
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.676e+00 2.280e+00 -1.174 0.242844
## casos.permil 2.956e-01 2.232e-02 13.244 < 2e-16 ***
## Pobla80 6.935e-02 8.350e-02 0.830 0.407998
## Pobla65 2.770e-02 2.760e-02 1.003 0.317740
## PoblaMid 2.612e-03 1.188e-02 0.220 0.826429
## PoblaData 2.796e-02 1.835e-02 1.524 0.130248
## PoblaDens -1.873e-02 4.769e-03 -3.927 0.000147 ***
## Mujeres 3.448e-02 1.597e-02 2.159 0.032948 *
## Urbano -5.229e-04 2.087e-03 -0.251 0.802598
## ExpectVida 5.198e-03 1.599e-02 0.325 0.745660
## NeontlMort -4.977e-03 7.136e-03 -0.697 0.486952
## DisMort 1.696e-03 1.172e-02 0.145 0.885178
## Lesion 4.963e-03 1.010e-02 0.491 0.624146
## EnfTrans 6.911e-03 7.071e-03 0.977 0.330411
## PBI 3.723e-03 2.544e-03 1.463 0.146102
## Tuberculosis -7.847e-05 3.133e-04 -0.250 0.802682
## Diabetes -4.676e-03 9.750e-03 -0.480 0.632407
## Medicos -3.004e-02 4.705e-02 -0.639 0.524401
## Camas -5.655e-02 2.012e-02 -2.810 0.005825 **
## ImmunSaramp -5.811e-04 2.837e-03 -0.205 0.838076
## TempMarzo -2.979e-03 4.790e-03 -0.622 0.535304
## HipTen -8.175e-03 4.681e-03 -1.746 0.083439 .
## BCG -3.779e-01 1.397e-01 -2.706 0.007848 **
## Tiempo 3.096e-03 1.903e-03 1.627 0.106491
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3067 on 115 degrees of freedom
## Multiple R-squared: 0.8632, Adjusted R-squared: 0.8358
## F-statistic: 31.55 on 23 and 115 DF, p-value: < 2.2e-16wise_modelo_muertes_permil_V4 = step(modelo_l10muertes_outliers, direction = c("both"),
trace = F)
formula_elegida = wise_modelo_muertes_permil_V4$call$formula
summary(wise_modelo_muertes_permil_V4)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaData +
## PoblaDens + Mujeres + EnfTrans + PBI + Camas + HipTen + BCG +
## Tiempo, data = df7[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7885 -0.1882 0.0072 0.1472 0.7662
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.162721 0.624221 -3.465 0.000724 ***
## casos.permil 0.294724 0.019165 15.378 < 2e-16 ***
## Pobla65 0.048674 0.008718 5.583 1.37e-07 ***
## PoblaData 0.024752 0.016943 1.461 0.146528
## PoblaDens -0.019352 0.004019 -4.815 4.10e-06 ***
## Mujeres 0.032079 0.011820 2.714 0.007571 **
## EnfTrans 0.004834 0.002032 2.379 0.018839 *
## PBI 0.004142 0.002167 1.911 0.058218 .
## Camas -0.052347 0.015767 -3.320 0.001175 **
## HipTen -0.009433 0.004136 -2.280 0.024251 *
## BCG -0.401135 0.125518 -3.196 0.001759 **
## Tiempo 0.002724 0.001681 1.621 0.107580
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2956 on 127 degrees of freedom
## Multiple R-squared: 0.8596, Adjusted R-squared: 0.8475
## F-statistic: 70.69 on 11 and 127 DF, p-value: < 2.2e-16Hemos logrado encontrar un mejor modelo, aumentando el r2-ajustado en 4.5% (0.8116 vs 0.8475) y, además, no hemos dropeado ninguna instancia, lo que es súmamente positivo. Vamos a ver los gráfico de residuos para asegurarnos de que no haya valores atipicos.
plot(wise_modelo_muertes_permil_V4)
ks.test(wise_modelo_muertes_permil_V4$residuals, "pnorm", mean(wise_modelo_muertes_permil_V4$residuals),
sd(wise_modelo_muertes_permil_V4$residuals))##
## One-sample Kolmogorov-Smirnov test
##
## data: wise_modelo_muertes_permil_V4$residuals
## D = 0.065046, p-value = 0.5988
## alternative hypothesis: two-sidedPareciera que tenemos distribucion normal y homocedasticidad en los residuos y no tenemos presencia de ningún outlier. Por último, vamos a correr un modelo robusto para comparar los beta con el modelo OLS y para ver el gráfico de residuos vs r-weigths
modelo_MM_covid_log = lmrobdetMM(formula_elegida,
data = select_if(df7,is.numeric), control = control)
summary(modelo_MM_covid_log)##
## Call:
## lmrobdetMM(formula = formula_elegida, data = select_if(df7, is.numeric),
## control = control)
## Residuals:
## Min 1Q Median 3Q Max
## -0.79458 -0.17096 0.01575 0.16144 0.85128
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.972730 0.690760 -2.856 0.00502 **
## casos.permil 0.275112 0.020656 13.319 < 2e-16 ***
## Pobla65 0.049847 0.009160 5.442 2.61e-07 ***
## PoblaData 0.025107 0.017233 1.457 0.14761
## PoblaDens -0.018085 0.004013 -4.506 1.48e-05 ***
## Mujeres 0.028493 0.013021 2.188 0.03048 *
## EnfTrans 0.005004 0.002146 2.331 0.02132 *
## PBI 0.004363 0.002297 1.900 0.05975 .
## Camas -0.050766 0.016571 -3.064 0.00267 **
## HipTen -0.008630 0.004282 -2.016 0.04594 *
## BCG -0.448764 0.133414 -3.364 0.00102 **
## Tiempo 0.003983 0.001782 2.235 0.02714 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.2881
## Multiple R-squared: 0.8604, Adjusted R-squared: 0.8483
## Convergence in 35 IRWLS iterationsplot(modelo_MM_covid_log$residuals,modelo_MM_covid_log$rweights,
main = "REGRESION ROBUSTA 3", xlab = "RESIDUOS",
ylab = "RWEIGHTS",col = "blue")
Claramente no se encuentra presencia de outliers y los beta estimados no presentan cambios significativos, por lo que podemos concluír en que nuestro modelo no se encuentra afectado por outliers significativos. Un punto más en cuanto a la selección del modelo, vemos que el método robusto mejora el r2-ajustado y el RSS vs la versión de OLS pero, como sabemos, este tipo de método tiene intervalos de confianza mayores, por lo que (al ser marginal la mejora), perferimos quedarnos con el modelo OLS.
final_model = wise_modelo_muertes_permil_V4
summary(lm.covid)##
## Call:
## lm(formula = l10muertes.permil ~ PoblaDens + Pobla80 + Urbano +
## Tuberculosis + Camas + TempMarzo + PBI, data = COVID)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08088 -0.29012 -0.03314 0.27929 1.29267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5981046 0.2012630 2.972 0.00352 **
## PoblaDens -0.0111873 0.0063386 -1.765 0.07990 .
## Pobla80 0.2040913 0.0338369 6.032 1.55e-08 ***
## Urbano 0.0082702 0.0025767 3.210 0.00167 **
## Tuberculosis -0.0005815 0.0003431 -1.695 0.09246 .
## Camas -0.1101291 0.0263087 -4.186 5.17e-05 ***
## TempMarzo -0.0143117 0.0057068 -2.508 0.01337 *
## PBI 0.0062466 0.0027181 2.298 0.02314 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4765 on 131 degrees of freedom
## Multiple R-squared: 0.6239, Adjusted R-squared: 0.6038
## F-statistic: 31.04 on 7 and 131 DF, p-value: < 2.2e-16summary(final_model)##
## Call:
## lm(formula = l10muertes.permil ~ casos.permil + Pobla65 + PoblaData +
## PoblaDens + Mujeres + EnfTrans + PBI + Camas + HipTen + BCG +
## Tiempo, data = df7[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7885 -0.1882 0.0072 0.1472 0.7662
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.162721 0.624221 -3.465 0.000724 ***
## casos.permil 0.294724 0.019165 15.378 < 2e-16 ***
## Pobla65 0.048674 0.008718 5.583 1.37e-07 ***
## PoblaData 0.024752 0.016943 1.461 0.146528
## PoblaDens -0.019352 0.004019 -4.815 4.10e-06 ***
## Mujeres 0.032079 0.011820 2.714 0.007571 **
## EnfTrans 0.004834 0.002032 2.379 0.018839 *
## PBI 0.004142 0.002167 1.911 0.058218 .
## Camas -0.052347 0.015767 -3.320 0.001175 **
## HipTen -0.009433 0.004136 -2.280 0.024251 *
## BCG -0.401135 0.125518 -3.196 0.001759 **
## Tiempo 0.002724 0.001681 1.621 0.107580
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2956 on 127 degrees of freedom
## Multiple R-squared: 0.8596, Adjusted R-squared: 0.8475
## F-statistic: 70.69 on 11 and 127 DF, p-value: < 2.2e-16plot(lm.covid)
plot(final_model)
Performance: vemos el salto significativo que hay entre los 2 modelos, logrando en el segundo un incremento en la performance marcado por: > un aumento del r2-ajustado de 0.6038 a 0.8475 (40%) > un descenso del RSS de 0.4765 a 0.2956 (62%) > la eliminación de la multicolinealidad > el manejo de outliers.
Cambios de variables: es destacable el cambio de variables seleccionadas y betas estimados, en el segundo modelo vemos que casos.permil es una variable preponderante (por ser la variable más grande y tener a su vez el beta estimado más grande) y en el primer modelo no la habíamos incluído. Algo parecido ocurre con la variable BCG, que toma reelevancia en el segundo modelo y en el primero no estaba.
par(mfrow=c(1, 2))
int_prediccion = predict(lm.covid, newdata = COVID, interval = c("predict"))
int_prediccion_2 = predict(wise_modelo_muertes_permil_V4, newdata = df7, interval = c("predict"))
plot(y = lm.covid$fitted.values,
x = COVID$l10muertes.permil,
col = "darkblue", ylim = c(-2, 5), cex = .8,
xlab = "VALORES OBSERVADOS", ylab = "VALORES PREDICHOS", main='Intervalos de predicción modelo 1')
# Intervalos de prediccion
points(x = COVID$l10muertes.permil,
y = int_prediccion[, 2], col = "darkred", lwd =3, cex = 0.5)
points(x = COVID$l10muertes.permil,
y = int_prediccion[, 3], col = "darkred", lwd =3, cex = 0.5)
legend("topleft", c("INTERVALO DE PREDICCION", "OBSERVACION"), fill = c("darkred", "darkblue"))
plot(y = wise_modelo_muertes_permil_V4$fitted.values,
x = df7$l10muertes.permil,
col = "darkblue", ylim = c(-2, 5), cex = .8,
xlab = "VALORES OBSERVADOS", ylab = "VALORES PREDICHOS", main='Intervalos de predicción modelo 2')
# Intervalos de prediccion
points(x = df7$l10muertes.permil,
y = int_prediccion_2[, 2], col = "darkred", lwd =3, cex = 0.5)
points(x = df7$l10muertes.permil,
y = int_prediccion_2[, 3], col = "darkred", lwd =3, cex = 0.5)
legend("topleft", c("INTERVALO DE PREDICCION", "OBSERVACION"), fill = c("darkred", "darkblue"))
En cuanto a los intervalos de predicción, vemos como en nuestro último modelo estos de vuelven más estrechos, acercándose a los valores reales de nuestra variable respuesta.
plot(final_model)
Observamos que un factor muy ponderante para disminuir las muertes debidas al COVID19, fue la estrategia de INMUNIZACION, por lo que sugeririamos a las autoridades que opten por una vacunacion TOTAL y no sectorizada, y se difundan mas publicidades para concientizar sobre los beneficios de vacunarse. Tambien aumentar la cantidad de campañas y centros de vacunacion; incluso de ser posible, dar incentivos a los que aun conservan dudas sobre si vacunarse, para que lo hagan.
Otro factor que influyo en gran medida fue la cantidad de casos por cada mil habitantes, por lo que recomendamos avanzar con estrategias de aislamiento preventivo y campañas de concientización de los cuidados escenciales como el uso de barbijo, higiene correcta de manos y distancia social.
Por último, tambien observamos que la poblacion mas afectada por el virus, fue la de mayor edad (mas de 65 años). Sugerimos implementar mas politicas de concientizacion a la sociedad sobre la vulneravilidad de este sector, para que se extremen más los recaudos. Tambien sugerimos agrandar los programas publicos de asistencia para personas mayores y lanzar voluntariados a la sociedad en general, para suplir esta demanda. De esta forma, la persona mayor puede disminuir la chance de contagio, evitando ir a lugares con mucha densidad de gente (supermercados, bancos, farmacias, etc), para cubrir sus necesidades basicas.